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Abstract. Numerical simulations are presented to study the stability of gaps opened by 
giant planets in 3D self-gravitating disks. In weakly self- gravitating disks, a few vortices 
develop at the gap edge and merge on orbital time- scales. The result is one large but weak 
vortex with Rossby number -0.01. In moderately self-gravitating disks, more vortices 
develop and their merging is resisted on dynamical time-scales. Self-gravity can sustain 
multi-vortex configurations, with Rossby number -0.2 to -0.1, over a time-scale of order 
100 orbits. Self-gravity also enhances the vortex vertical density stratification, even in 
disks with initial Toomre parameter of order 10. However, vortex formation is suppressed 
in strongly self-gravitating disks and replaced by a global spiral instability associated with 
the gap edge which develops during gap formation. 



1 Introduction 



Gaps induced by planets in protoplanetary disks can become dynamically unstable if the disk viscosity 

is sufficiently small | IP. This is because planetary gap edges are associated with potential vorticity 

. . or vortensity extrema Q2L IsL the existence of which is necessary for instability [4]. Gap edges may 

^ ^ undergo vortex formation in weakly self-gravitating disks (associated with vortensity minima) or a 

k> ' spiral instability (associated with vortensity maxima) in strongly self-gravitating yet Toomre- stable 

$_ j .disks 13, |a]. Development of such instabilities can significantly affect planetary migration [7] and 

Cy _"dust evolution QSQ - These studies have employed 2D disk models, but gap edges have characteristic 

widths of the disk local scale-height, so it is necessary to extend the study of gap stability to 3D. 

2 Numerical simulations with ZEUS-MP 

The system is an inviscid, non-magnetized 3D fluid disk embedded with a giant planet of mass M p , 
both rotating about a central star of mass M*. Spherical co-ordinates (r, 6, (p) centered about the star 
are adopted. Units are such that G = M* = 1 , where G is the gravitational constant. 

The disk is governing by the Euler equations coupled with self-gravity through the Poisson equa- 
tion. The equation of state is locally isothermal, so the sound-speed is c s = H£lk, where H = hR is the 
isothermal scale-height with constant aspect-ratio h, £l 2 k = GM*/R 3 and R = rsinO. Each disk model 
is labelled by its minimum Keplerian Toomre parameter Qq, located at the outer disk boundary. The 
planet is regarded as an external potential and held on a circular Keplerian orbit of radius r p at the 
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Figure 1. Relative density perturbation in the (r, 6) plane, chosen at the azimuth that intercepts the vortex centroid 
formed in two disk models with go = 8. The perturbed meridional flow is also shown. 



midplane. Time is quoted in units of Po = 2^/Q^(r /7 ). The Hill radius Vh = (M p /3M*) l/3 r p is used in 
some of the plots. 

The self-gravitating hydrodynamic equations are evolved with the ZEUS -MP finite difference code 
I9D. The computational domain, unless otherwise stated, is r e [1, 25], € [#min, 7*72] and (p € [0, 2n] 
where tan (n/2 - 0^) = 2h. Boundary conditions are outflow in r, reflecting in 6 and periodic in 
0. The numerical resolution is N r x Ng x N$ = 256 x 32 x 512. See 110] for further details for the 
simulation setup. 



3 Results 

3.1 Weakly self-gravitating disks (g = 8) 

Two simulations with go - 8 were run, one with self-gravity and the other without. M p = 0.002M* 
and h = 0.07 are adopted. Both cases developed 2 — 3 vortices early on but the quasi- steady state is a 

single vortex with Ro 0.01, where the Rossby number is defined as Ro = a> z /(2£l), where co z is 

the absolute vertical vorticity and (CI) is the azimuthally-averaged angular speed. 

The vortices differ noticeably in the (r, 6) plane. This is shown in Fig. [T] Self-gravity enhances the 
vertical density stratification of the vortex, with the midplane density enhancement being ^ 50% larger 
than that in the non- self-gravitating run. The initial Keplerian Toomre parameter at the radius of vortex 
formation is about 10, but even this is sufficient to affect the vortex vertical structure. The creation of 
vortensity minima lowers the Toomre parameter, which further decreases with vortex formation since 
they are over-densities. Thus, self-gravity can become important in the perturbed state even if it is 
negligible initially. 



3.2 Moderately self-gravitating disks (go = 3) 

Fig. [2] shows the relative density perturbation and Rossby number at the end of the a simulation with 
Qo = 3. Self-gravity is included. The 5-vortex configuration is sustained from its initial development, 
unlike in the weakly self- gravitating case where merging occurred over the same time- scale. The 
preference for linear vortex modes with higher azimuthal wavenumber m with increasing strength of 
self-gravity was observed in 2D simulations iQ Hill , and persists in 3D. The smaller vortices here 

are stronger than the single vortex in previous case, with Rossby number Ro 0.2 and the relative 

density perturbation has significant vertical dependence. 
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Figure 2. Multi-vortex configuration at the end of a simulation for Q = 3 (t = 50P X with h = 0.07 and 
M p - 0.002M* . The Rossby number and relative density perturbation in the (r, 0) plane at 6 - n/2 are shown on 
the left; and in the (r, 6) plane at = O on the right, where O is the vortex azimuth denoted by dotted lines. 




Figure 3. Simulation with Q = 3.0, h = 0.05 and M p = 10 3 M*. The midplane relative density perturbation is 
shown at t - 300P (left), t - 405P (middle) and t - 505P (right). The minimum Rossby number was found 
to be Ro = -0.11 (left), Ro = -0.09 (middle) and Ro = +0.03 (right). 



3.2. 1 Long term simulation 

A smaller disk model, with r e [2, 20], was simulated up to t ~ 500Po- M p = 10" 3 M* and h = 0.05 
were adopted for this run. Fig. shows the relative density perturbation towards the end of the 
simulation. The multi- vortex configuration lasted ~ 200 orbits at the vortex radius. Notice a vortex 
may reach comparable over-densities to the final post-merger vortex in the weakly self-gravitating 
disk. It was observed that \Ro\ decreased from 0.2 at the onset of vortex formation, to 0.1 towards the 
end of the simulation, which may be due to limited numerical resolution. 



3.3 Gap edge spiral instability (2 = 1.5) 

The linear vortex instability can be suppressed by strong self-gravity. To demonstrate this, a disk 
model with go = 1.5, h = 0.05 and M p = 10" 3 M* was simulated. Fig.|4]shows the development of an 
m = 2 spiral mode associated with the outer gap edge. This instability occurs during gap formation 
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Figure 4. Spiral instability associated with the outer gap edge opened by a giant planet (not a classic Toomre 



instability). The disk model is go 
suppress vortex formation. 



1.5, h - 0.05 with M p = 10 M*. Self-gravity is sufficiently strong to 



and supplies positive disk-on-planet co-orbital torques, because the over-density protrudes the outer 
gap edge and approaches the planet from upstream. The disturbance is significantly stratified, with 
most of the perturbation confined near the midplane. The global spiral pattern appears transient, 
having decreased in amplitude by t = 50Pq, but this is likely a radial boundary condition effect. 

4 Discussion 

Direct numerical simulations of 3D self- gravitatin g d isk-planet systems confirm the stability proper- 
ties of gap edges previously explored in 2D [^ISUJIlLJ]- Vertical self-gravity enhances the density 
stratification of a vortex. Given the vortex instability is only expected to occur in low viscosity regions 
of protoplanetary disks — dead zones — which are overlaid by actively accreting layers [J13I1 . it may 
be advantageous to have self-gravity confining the over-density near the midplane, thereby mitigate 
upper disk boundary effects and make the instability a more robust mechanism for vortex formation. 
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